home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 1 Issue 2 / PDCD-1 - Issue 02.iso / _utilities / utilities / 001 / meschach / !Meschach / c / zschur < prev    next >
Text File  |  1994-01-13  |  11KB  |  376 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26. /*    
  27.     File containing routines for computing the Schur decomposition
  28.     of a complex non-symmetric matrix
  29.     See also: hessen.c
  30.     Complex version
  31. */
  32.  
  33.  
  34. #include    <stdio.h>
  35. #include    <math.h>
  36. #include    "zmatrix.h"
  37. #include        "zmatrix2.h"
  38.  
  39.  
  40. #define    is_zero(z)    ((z).re == 0.0 && (z).im == 0.0)
  41. #define    b2s(t_or_f)    ((t_or_f) ? "TRUE" : "FALSE")
  42.  
  43.  
  44. /* zschur -- computes the Schur decomposition of the matrix A in situ
  45.     -- optionally, gives Q matrix such that Q^T.A.Q is upper triangular
  46.     -- returns upper triangular Schur matrix */
  47. ZMAT    *zschur(A,Q)
  48. ZMAT    *A, *Q;
  49. {
  50.     int        i, j, iter, k, k_min, k_max, k_tmp, n, split;
  51.     Real    c;
  52.     complex    det, discrim, lambda, lambda0, lambda1, s, sum, ztmp;
  53.     complex    x, y;    /* for chasing algorithm */
  54.     complex    **A_me;
  55.     static    ZVEC    *diag=ZVNULL;
  56.     
  57.     if ( ! A )
  58.     error(E_NULL,"zschur");
  59.     if ( A->m != A->n || ( Q && Q->m != Q->n ) )
  60.     error(E_SQUARE,"zschur");
  61.     if ( Q != ZMNULL && Q->m != A->m )
  62.     error(E_SIZES,"zschur");
  63.     n = A->n;
  64.     diag = zv_resize(diag,A->n);
  65.     MEM_STAT_REG(diag,TYPE_ZVEC);
  66.     /* compute Hessenberg form */
  67.     zHfactor(A,diag);
  68.     
  69.     /* save Q if necessary, and make A explicitly Hessenberg */
  70.     zHQunpack(A,diag,Q,A);
  71.  
  72.     k_min = 0;    A_me = A->me;
  73.  
  74.     while ( k_min < n )
  75.     {
  76.     /* find k_max to suit:
  77.        submatrix k_min..k_max should be irreducible */
  78.     k_max = n-1;
  79.     for ( k = k_min; k < k_max; k++ )
  80.         if ( is_zero(A_me[k+1][k]) )
  81.         {    k_max = k;    break;    }
  82.  
  83.     if ( k_max <= k_min )
  84.     {
  85.         k_min = k_max + 1;
  86.         continue;        /* outer loop */
  87.     }
  88.  
  89.     /* now have r x r block with r >= 2:
  90.        apply Francis QR step until block splits */
  91.     split = FALSE;        iter = 0;
  92.     while ( ! split )
  93.     {
  94.         complex    a00, a01, a10, a11;
  95.         iter++;
  96.         
  97.         /* set up Wilkinson/Francis complex shift */
  98.         /* use the smallest eigenvalue of the bottom 2 x 2 submatrix */
  99.         k_tmp = k_max - 1;
  100.  
  101.         a00 = A_me[k_tmp][k_tmp];
  102.         a01 = A_me[k_tmp][k_max];
  103.         a10 = A_me[k_max][k_tmp];
  104.         a11 = A_me[k_max][k_max];
  105.         ztmp.re = 0.5*(a00.re - a11.re);
  106.         ztmp.im = 0.5*(a00.im - a11.im);
  107.         discrim = zsqrt(zadd(zmlt(ztmp,ztmp),zmlt(a01,a10)));
  108.         sum.re  = 0.5*(a00.re + a11.re);
  109.         sum.im  = 0.5*(a00.im + a11.im);
  110.         lambda0 = zadd(sum,discrim);
  111.         lambda1 = zsub(sum,discrim);
  112.         det = zsub(zmlt(a00,a11),zmlt(a01,a10));
  113.         if ( zabs(lambda0) > zabs(lambda1) )
  114.         lambda = zdiv(det,lambda0);
  115.         else
  116.         lambda = zdiv(det,lambda1);
  117.  
  118.         /* perturb shift if convergence is slow */
  119.         if ( (iter % 10) == 0 )
  120.         {
  121.         lambda.re += iter*0.02;
  122.         lambda.im += iter*0.02;
  123.         }
  124.  
  125.         /* set up Householder transformations */
  126.         k_tmp = k_min + 1;
  127.  
  128.         x = zsub(A->me[k_min][k_min],lambda);
  129.         y = A->me[k_min+1][k_min];
  130.  
  131.         /* use Givens' rotations to "chase" off-Hessenberg entry */
  132.         for ( k = k_min; k <= k_max-1; k++ )
  133.         {
  134.         zgivens(x,y,&c,&s);
  135.         zrot_cols(A,k,k+1,c,s,A);
  136.         zrot_rows(A,k,k+1,c,s,A);
  137.         if ( Q != ZMNULL )
  138.             zrot_cols(Q,k,k+1,c,s,Q);
  139.  
  140.         /* zero things that should be zero */
  141.         if ( k > k_min )
  142.             A->me[k+1][k-1].re = A->me[k+1][k-1].im = 0.0;
  143.  
  144.         /* get next entry to chase along sub-diagonal */
  145.         x = A->me[k+1][k];
  146.         if ( k <= k_max - 2 )
  147.             y = A->me[k+2][k];
  148.         else
  149.             y.re = y.im = 0.0;
  150.         }
  151.  
  152.         for ( k = k_min; k <= k_max-2; k++ )
  153.         {
  154.         /* zero appropriate sub-diagonals */
  155.         A->me[k+2][k].re = A->me[k+2][k].im = 0.0;
  156.         }
  157.  
  158.         /* test to see if matrix should split */
  159.         for ( k = k_min; k < k_max; k++ )
  160.         if ( zabs(A_me[k+1][k]) < MACHEPS*
  161.             (zabs(A_me[k][k])+zabs(A_me[k+1][k+1])) )
  162.         {
  163.             A_me[k+1][k].re = A_me[k+1][k].im = 0.0;
  164.             split = TRUE;
  165.         }
  166.  
  167.     }
  168.     }
  169.     
  170.     /* polish up A by zeroing strictly lower triangular elements
  171.        and small sub-diagonal elements */
  172.     for ( i = 0; i < A->m; i++ )
  173.     for ( j = 0; j < i-1; j++ )
  174.         A_me[i][j].re = A_me[i][j].im = 0.0;
  175.     for ( i = 0; i < A->m - 1; i++ )
  176.     if ( zabs(A_me[i+1][i]) < MACHEPS*
  177.         (zabs(A_me[i][i])+zabs(A_me[i+1][i+1])) )
  178.         A_me[i+1][i].re = A_me[i+1][i].im = 0.0;
  179.  
  180.     return A;
  181. }
  182.  
  183.  
  184. #if 0
  185. /* schur_vecs -- returns eigenvectors computed from the real Schur
  186.         decomposition of a matrix
  187.     -- T is the block upper triangular Schur matrix
  188.     -- Q is the orthognal matrix where A = Q.T.Q^T
  189.     -- if Q is null, the eigenvectors of T are returned
  190.     -- X_re is the real part of the matrix of eigenvectors,
  191.         and X_im is the imaginary part of the matrix.
  192.     -- X_re is returned */
  193. MAT    *schur_vecs(T,Q,X_re,X_im)
  194. MAT    *T, *Q, *X_re, *X_im;
  195. {
  196.     int    i, j, limit;
  197.     Real    t11_re, t11_im, t12, t21, t22_re, t22_im;
  198.     Real    l_re, l_im, det_re, det_im, invdet_re, invdet_im,
  199.         val1_re, val1_im, val2_re, val2_im,
  200.         tmp_val1_re, tmp_val1_im, tmp_val2_re, tmp_val2_im, **T_me;
  201.     Real    sum, diff, discrim, magdet, norm, scale;
  202.     static VEC    *tmp1_re=VNULL, *tmp1_im=VNULL,
  203.             *tmp2_re=VNULL, *tmp2_im=VNULL;
  204.  
  205.     if ( ! T || ! X_re )
  206.         error(E_NULL,"schur_vecs");
  207.     if ( T->m != T->n || X_re->m != X_re->n ||
  208.         ( Q != MNULL && Q->m != Q->n ) ||
  209.         ( X_im != MNULL && X_im->m != X_im->n ) )
  210.         error(E_SQUARE,"schur_vecs");
  211.     if ( T->m != X_re->m ||
  212.         ( Q != MNULL && T->m != Q->m ) ||
  213.         ( X_im != MNULL && T->m != X_im->m ) )
  214.         error(E_SIZES,"schur_vecs");
  215.  
  216.     tmp1_re = v_resize(tmp1_re,T->m);
  217.     tmp1_im = v_resize(tmp1_im,T->m);
  218.     tmp2_re = v_resize(tmp2_re,T->m);
  219.     tmp2_im = v_resize(tmp2_im,T->m);
  220.     MEM_STAT_REG(tmp1_re,TYPE_VEC);
  221.     MEM_STAT_REG(tmp1_im,TYPE_VEC);
  222.     MEM_STAT_REG(tmp2_re,TYPE_VEC);
  223.     MEM_STAT_REG(tmp2_im,TYPE_VEC);
  224.  
  225.     T_me = T->me;
  226.     i = 0;
  227.     while ( i < T->m )
  228.     {
  229.         if ( i+1 < T->m && T->me[i+1][i] != 0.0 )
  230.         {    /* complex eigenvalue */
  231.         sum  = 0.5*(T_me[i][i]+T_me[i+1][i+1]);
  232.         diff = 0.5*(T_me[i][i]-T_me[i+1][i+1]);
  233.         discrim = diff*diff + T_me[i][i+1]*T_me[i+1][i];
  234.         l_re = l_im = 0.0;
  235.         if ( discrim < 0.0 )
  236.         {    /* yes -- complex e-vals */
  237.             l_re = sum;
  238.             l_im = sqrt(-discrim);
  239.         }
  240.         else /* not correct Real Schur form */
  241.             error(E_RANGE,"schur_vecs");
  242.         }
  243.         else
  244.         {
  245.         l_re = T_me[i][i];
  246.         l_im = 0.0;
  247.         }
  248.  
  249.         v_zero(tmp1_im);
  250.         v_rand(tmp1_re);
  251.         sv_mlt(MACHEPS,tmp1_re,tmp1_re);
  252.  
  253.         /* solve (T-l.I)x = tmp1 */
  254.         limit = ( l_im != 0.0 ) ? i+1 : i;
  255.         /* printf("limit = %d\n",limit); */
  256.         for ( j = limit+1; j < T->m; j++ )
  257.         tmp1_re->ve[j] = 0.0;
  258.         j = limit;
  259.         while ( j >= 0 )
  260.         {
  261.         if ( j > 0 && T->me[j][j-1] != 0.0 )
  262.         {   /* 2 x 2 diagonal block */
  263.             /* printf("checkpoint A\n"); */
  264.             val1_re = tmp1_re->ve[j-1] -
  265.               __ip__(&(tmp1_re->ve[j+1]),&(T->me[j-1][j+1]),limit-j);
  266.             /* printf("checkpoint B\n"); */
  267.             val1_im = tmp1_im->ve[j-1] -
  268.               __ip__(&(tmp1_im->ve[j+1]),&(T->me[j-1][j+1]),limit-j);
  269.             /* printf("checkpoint C\n"); */
  270.             val2_re = tmp1_re->ve[j] -
  271.               __ip__(&(tmp1_re->ve[j+1]),&(T->me[j][j+1]),limit-j);
  272.             /* printf("checkpoint D\n"); */
  273.             val2_im = tmp1_im->ve[j] -
  274.               __ip__(&(tmp1_im->ve[j+1]),&(T->me[j][j+1]),limit-j);
  275.             /* printf("checkpoint E\n"); */
  276.             
  277.             t11_re = T_me[j-1][j-1] - l_re;
  278.             t11_im = - l_im;
  279.             t22_re = T_me[j][j] - l_re;
  280.             t22_im = - l_im;
  281.             t12 = T_me[j-1][j];
  282.             t21 = T_me[j][j-1];
  283.  
  284.             scale =  fabs(T_me[j-1][j-1]) + fabs(T_me[j][j]) +
  285.             fabs(t12) + fabs(t21) + fabs(l_re) + fabs(l_im);
  286.  
  287.             det_re = t11_re*t22_re - t11_im*t22_im - t12*t21;
  288.             det_im = t11_re*t22_im + t11_im*t22_re;
  289.             magdet = det_re*det_re+det_im*det_im;
  290.             if ( sqrt(magdet) < MACHEPS*scale )
  291.             {
  292.                 det_re = MACHEPS*scale;
  293.             magdet = det_re*det_re+det_im*det_im;
  294.             }
  295.             invdet_re =   det_re/magdet;
  296.             invdet_im = - det_im/magdet;
  297.             tmp_val1_re = t22_re*val1_re-t22_im*val1_im-t12*val2_re;
  298.             tmp_val1_im = t22_im*val1_re+t22_re*val1_im-t12*val2_im;
  299.             tmp_val2_re = t11_re*val2_re-t11_im*val2_im-t21*val1_re;
  300.             tmp_val2_im = t11_im*val2_re+t11_re*val2_im-t21*val1_im;
  301.             tmp1_re->ve[j-1] = invdet_re*tmp_val1_re -
  302.                     invdet_im*tmp_val1_im;
  303.             tmp1_im->ve[j-1] = invdet_im*tmp_val1_re +
  304.                     invdet_re*tmp_val1_im;
  305.             tmp1_re->ve[j]   = invdet_re*tmp_val2_re -
  306.                     invdet_im*tmp_val2_im;
  307.             tmp1_im->ve[j]   = invdet_im*tmp_val2_re +
  308.                     invdet_re*tmp_val2_im;
  309.             j -= 2;
  310.             }
  311.             else
  312.         {
  313.             t11_re = T_me[j][j] - l_re;
  314.             t11_im = - l_im;
  315.             magdet = t11_re*t11_re + t11_im*t11_im;
  316.             scale = fabs(T_me[j][j]) + fabs(l_re);
  317.             if ( sqrt(magdet) < MACHEPS*scale )
  318.             {
  319.                 t11_re = MACHEPS*scale;
  320.             magdet = t11_re*t11_re + t11_im*t11_im;
  321.             }
  322.             invdet_re =   t11_re/magdet;
  323.             invdet_im = - t11_im/magdet;
  324.             /* printf("checkpoint F\n"); */
  325.             val1_re = tmp1_re->ve[j] -
  326.               __ip__(&(tmp1_re->ve[j+1]),&(T->me[j][j+1]),limit-j);
  327.             /* printf("checkpoint G\n"); */
  328.             val1_im = tmp1_im->ve[j] -
  329.               __ip__(&(tmp1_im->ve[j+1]),&(T->me[j][j+1]),limit-j);
  330.             /* printf("checkpoint H\n"); */
  331.             tmp1_re->ve[j] = invdet_re*val1_re - invdet_im*val1_im;
  332.             tmp1_im->ve[j] = invdet_im*val1_re + invdet_re*val1_im;
  333.             j -= 1;
  334.         }
  335.         }
  336.  
  337.         norm = v_norm_inf(tmp1_re) + v_norm_inf(tmp1_im);
  338.         sv_mlt(1/norm,tmp1_re,tmp1_re);
  339.         if ( l_im != 0.0 )
  340.         sv_mlt(1/norm,tmp1_im,tmp1_im);
  341.         mv_mlt(Q,tmp1_re,tmp2_re);
  342.         if ( l_im != 0.0 )
  343.         mv_mlt(Q,tmp1_im,tmp2_im);
  344.         if ( l_im != 0.0 )
  345.         norm = sqrt(in_prod(tmp2_re,tmp2_re)+in_prod(tmp2_im,tmp2_im));
  346.         else
  347.         norm = v_norm2(tmp2_re);
  348.         sv_mlt(1/norm,tmp2_re,tmp2_re);
  349.         if ( l_im != 0.0 )
  350.         sv_mlt(1/norm,tmp2_im,tmp2_im);
  351.  
  352.         if ( l_im != 0.0 )
  353.         {
  354.         if ( ! X_im )
  355.         error(E_NULL,"schur_vecs");
  356.         set_col(X_re,i,tmp2_re);
  357.         set_col(X_im,i,tmp2_im);
  358.         sv_mlt(-1.0,tmp2_im,tmp2_im);
  359.         set_col(X_re,i+1,tmp2_re);
  360.         set_col(X_im,i+1,tmp2_im);
  361.         i += 2;
  362.         }
  363.         else
  364.         {
  365.         set_col(X_re,i,tmp2_re);
  366.         if ( X_im != MNULL )
  367.             set_col(X_im,i,tmp1_im);    /* zero vector */
  368.         i += 1;
  369.         }
  370.     }
  371.  
  372.     return X_re;
  373. }
  374.  
  375. #endif
  376.